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We report quantum Monte Carlo calculations of biexciton binding energies in ideal two- 
dimensional bilayer systems with isotropic electron and hole masses. We have also calculated 
exciton-exciton interaction potentials, and pair distribution functions for electrons and holes in 
bound biexcitons. Comparing our data with results obtained in a recent study using a model 
exciton-exciton potential [C. Schindler and R. Zimmermann, Phys. Rev. B 78, 045313 (2008)], we 
find a somewhat larger range of layer separations at which biexcitons are stable. We find that 
individual excitons retain their identity in bound biexcitons for large layer separations. 
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I. INTRODUCTION 

Electrons and holes in semiconduetors can combine to 
form hydrogen-like bound states called excitons. The cre- 
ation and recombination of excitons is one of the prin- 
cipal mechanisms by which light interacts with semicon- 
ductors. Furthermore, excitonic systems possess a num- 
ber of unusual properties, such as the ability to trans- 
port energy without transporting charge, suggesting a 
range of applications in novel electronic devices. Exci- 
tonic systems have therefore been the subject of numer- 
ous experimentalii^'^'^'^ and theoretical^ii'^ studies in re- 
cent years. 

In the low-density limit, excitons may be regarded as 
weakly-interacting neutral bosons Bose-Einstein con- 
densation of excitons is therefore possible. However, 
it has proved to be difficult to obtain exciton lifetimes 
that are sufficiently long for thermalization to take place, 
and the current experimental evidence for Bose-Einstein 
condensation is inconclusive.^^ One promising approach 
for overcoming the problem of short lifetimes is the de- 
velopment of coupled-quantum-well (bilayer) systems, in 
which thin layers of semiconductor and an applied elec- 
tric field in the growth direction are used to confine the 
electrons and holes to spatially separated, parallel, quasi- 
two-dimensional wells, hindering recombination and ex- 
tending exciton lifetimes fii^i^ 

At present our understanding of the exciton-exciton in- 
teraction in bilayer systems is limited. On the one hand 
there is a repulsive electrostatic interaction between ex- 
citons. For example, if the layer separation is nonzero 
then the excitons have parallel dipole moments, giving an 
asymptotically dominant repulsive interaction. Further- 
more, the static charge distribution of each exciton has a 
permanent quadrupole moment in general (even at zero 
layer separation, provided the electron and hole masses 
differ), giving another repulsive interaction term J. On the 
other hand, fluctuating dipole (van der Waals) forces re- 
sult in an attraction between excitons at short range. 
Because of the existence of the van der Waals forces, it is 
sometimes possible for biexcitons (bound states of pairs 
of excitons) to form. A better understanding of the in- 



teraction between excitons in coupled quantum wells will 
facilitate the interpretation of experimental data, in par- 
ticular enabling the determination of the exciton densi- 
ties achieved in experiments. 

The dependence of exciton and biexciton binding en- 
ergies on the layer separation has been investigated by 
Tan et ai, who found that, while the exciton binding en- 
ergy decays slowly as the inverse of the layer separation, 
the biexciton binding energy decays extremely rapidlyi^ 
Recent studies of the exciton-exciton interaction using a 
heavy-hole approximation have found there to be a criti- 
cal layer separation for each electron/hole mass ratio, be- 
yond which biexcitons become unstable with respect to 
dissociation into two separate excitons^iii In this article 
we report quantum Monte Carlo (QMC) calculations of 
the binding energies of biexcitons in bilayer systems and 
exciton-exciton interaction potentials. 

The rest of this article is arranged as follows. In Sec. 
im we describe our calculations of the binding energies 
of biexcitons and investigate the range of layer separa- 
tions and mass ratios for which biexcitons are stable. In 
Sec, mil we present our data for the exciton-exciton inter- 
action potential. In Sec. IIVI we report pair-distribution 
functions (PDF's) for biexcitons. Finally, we draw our 
conclusions in Sec. El We use Hartree atomic units 
(h = \e\ = rrie = 47reo = 1) throughout this article, 
although we report final energies in exciton Rydbergs 
{R* = /xe^/[2(47reoe)^?i^], where fi = rriemh/ {rrie + mfi) is 
the reduced mass of an exciton and rrie and m/i are the 
electron and hole masses) and lengths in terms of exciton 
Bohr radu [a*^ = 47reoe/i^/(/ie^)]. 



II. BIEXCITON BINDING ENERGIES 

We have modeled the coupled-quantum- well system by 
an idealized two-dimensional (2D) bilayer, in which the 
electrons and holes are confined to two parallel planes, 
and the effective mass tensors of the electrons and holes 
are isotropic. In reality, electrons and holes are free to 
move within quantum wells that are of finite width (e.g., 
one experimental setup*^ has well widths of 8 nm and a 
well separation of 4 nm), although the Coulomb attrac- 
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tion between electrons and holes should keep the particles 
confined to the inner edges of their respective wells. We 
have also restricted our attention to biexciton systems in 
which the two electrons have opposite spins, as do the 
two holes, because this is the ground-state spin configu- 
ration. 

We have studied biexcitons and exciton-exciton in- 
teractions using the variational and diffusion quantum 
Monte Carlo (VMC and DMC) methods. In the VMC 
method the expectation value of the Hamiltonian with 
respect to a trial wave function is calculated using a 
stochastic integration technique.— Trial wave functions 
usually contain a number of free parameters, to be opti- 
mized by minimizing either the energy expectation value 
or the variance of the energy. DMC is a stochastic pro- 
jector technique for solving the many-body Schrodinger 
equationj ^^i^^ DMC is in principle exact for systems with 
nodeless ground-state wave functions, such as the biexci- 
ton systems studied in this work. 

Our trial wave function was similar to that of Tan et 
al.^ with additional flexibility provided by multiplication 
by a two-body Jastrow factor This wave function sat- 
isfies the Kato cusp conditions when particles coincidcji^ 
and reduces to the form appropriate for two isolated exci- 
tons when the excitons are far apart. We also carried out 
some calculations using a three-body Jastrow factorJ^ 
For a typical case where i?vMC ~ E^mc = 3 x 10"'' i?*, 
the reduction in the VMC energy from the inclusion of a 
three-body term was 10~^i?,*. Obtaining the best possi- 
ble trial wave function was especially important for the 
PDF calculations described in Sec. IIVI 

We optimized the free parameters in our wave function 
by unreweighted variance minimizatio n^^i^^'^^ and linear- 
least-squares energy minimizatioui^ The trial wave func- 
tion can describe the dissociated system more accurately 
than it can describe the bound system; hence energy 
minimization is the more sensible choice for investigat- 
ing binding, although this depends upon initial param- 
eters and configurations. Although the DMC energy is 
independent of the trial wave function, the statistical effi- 
ciency of the method is increased when the wave function 
is improved. By using a more flexible wave function, we 
have been able to achieve considerably smaller error bars 
than Tan et alJ- All the QMC calculations reported in 
this article were performed using the CASINO programii^ 

For each layer separation d and electron/hole mass ra- 
tio cr, the biexciton binding energy was calculated as 
Eh ~ 2Ex — Exx, where Ex is the energy of a single 
exciton and Exx is the energy of the four-body biexci- 
ton system. The exciton energy Ex was obtained using 
a numerically exact Runge-Kutta integration technique 
as described in Ref. 0, while DMC was used to calculate 
the biexciton energy Exx- The DMC energies were con- 
verged with respect to time step and population size; any 
remaining bias is much smaller than the statistical error 
bars. 

Biexciton binding energies for a = 0.3, 0.5, and 1 are 
shown in Figs. [1] ^ and [31 respectively. It can be seen 



that our results are close to those of Tan et al.'r- the differ- 
ence arises from our use of exact single-exciton energies. 
Tan et al. used Eq. (3) of Ref. to generate Ex values, in- 
troducing a small, systematic error. Removing this error 
reveals that our DMC data are in statistical agreement 
with those of Tan et al. The random errors in our data 
are much smaller, so we can locate the layer separation 
at which the biexciton ceases to be bound. Tan et al. 
fltted an exponential form to their binding energy data, 
which resulted in the erroneous conclusion that biexciton 
binding persists to infinite layer separation. 

Our DMC results show some deviation from the bind- 
ing energies obtained by Schindler and Zimmermann;^ 
especially when c? — > 0, when E^ — > 0, and when the 
mass ratio is close to 1, because we have performed a full 
simulation of all four particles in the biexciton, whereas 
they simulated a pair of excitons interacting via a model 
potential. The deviation of our binding energies from 
those of Schindler and Zimmermann is approximately 
4 X 10"'^ i?* where E^ 0. At smaller d the agreement 
is much better, but below d « O.la^ we find larger dif- 
ferences, reaching a maximum of almost O.li?* at d = 0, 
as shown in Figs. [T] and [2] 

As can be seen in Fig. HI which shows the range of 
a and d over which the biexciton is stable, we find a 
somewhat larger region of stability for the biexciton than 
Schindler and Zimmermann. Let dcritlc) be the criti- 
cal layer separation, beyond which the biexciton is un- 
bound. As (T ^ 0, the heavy-hole approximation made 
by Schindler and Zimmermann becomes increasingly ac- 
curate, and our results for (icrit(O) agree with theirs. On 
the other hand, for a ~ I their interaction potential is 
less accurate and our value of dcrit(l) is therefore signif- 
icantly higher than theirs. 

Our data are mostly in excellent agreement with those 
of Meyertholen and Fogler, although at small cr we find a 
slightly larger region of biexciton stability. This is not an 
artifact of the extrapolation, which followed the scheme 
set out in Ref. [2ll, for we were able to find points with 
nonzero binding energies outside the region of stability 
defined by Meyertholen and Fogler. This is consistent 
with the variational principle that applies to their results. 

One may parameterize the boundary of the region of 
biexciton stability in Fig. 3) Expressing dcrit in terms 
of cr -I- cr~^ ensures that the correct behavior is ob- 
served upon exchanging the electron and hole masses [i.e., 
'^crit(cr^^) = rfcrit(cr)]- A Suitable fitting function is 



tanh 



G\/7 



+ 0.93, (1) 



where the parameter values F ~ 1.19(5) and G = 
—0.50(4) give a error of 0.4 per data point. The func- 
tional form of Eq. ([1]) satisfies most of the conditions 
derived in Ref. [2l|: d^j.;^(0) is infinite, d[.j.;j(l) = and 
rfcrit(O) - dcriticr) (X y/a foT a <^ 1. 
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FIG. 1: (Color online) Biexciton binding energy Ei, as a func- 
tion of layer separation d for electron/hole mass ratio a — 0.3. 
The upper panel shows the binding energy for layer separa- 
tions close to the critical separation; the lower panel shows 
the binding energy for small layer separations. 
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FIG. 2: (Color online) Biexciton binding energy as a func- 
tion of layer separation d for electron/hole mass ratio a = 0.5. 
The upper panel shows the binding energy for layer separa- 
tions close to the critical separation; the lower panel shows 
the binding energy for small layer separations. 



III. EXCITON-EXCITON INTERACTION 

The exciton-exciton interaction potential Ei(R) at 
separation R is defined to be the energy of a biexciton 
system in which the centers of mass of the two excitons 
are constrained to be a distance R apart, minus the en- 
ergies of two isolated excitons. The Hamiltonian for the 
constrained biexciton system may be written as 

2/i n r2 

1 1 

+ |R+;f;(-r2 + ri)| + |R+;^(-ri+r2)| 
1 1 

(2) 

where ri and r2 are the electron-hole separations within 
the two excitons. The first two potential terms repre- 
sent the intra-exciton electron-hole potentials, followed 
by the hole-hole, electron-electron, and finally the two 
inter-exciton electron-hole terms. DMC calculations can 



15 


' 1 ' 1 






* V 


• DMC (present work) 
Schindler & Zimmermann 

♦ DMC (Tan a/.) 






* V 






m 

'o 





















1 












0.26 0.28 


0.3 0.32 0.34 0.36 0.38 





FIG. 3: (Color online) Biexciton binding energy Eb as a func- 
tion of layer separation d for equal electron and hole masses 
{a — 1). The square shows Schindler and Zimmermann's es- 
timate of the critical point at which the biexciton ceases to 
be bound.— 
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FIG. 4: (Color online) The region of biexciton stability from 
DMC calculations compared with that found by Schindler 
and Zimmermann^ and Meyertholen and Fogler.— The criti- 
cal points were found by extrapolating the biexciton binding 
energies to zero using the fitting form set out in Ref. [2l|. The 
statistical errors are comparable to the size of the symbols. 



then be performed for an effective two-particle system, 
with coordinates ri and r2. The kinetic-energy operator 
only includes derivatives with respect to in-plane coordi- 
nates. The form of trial wave function was the same as 
that used in Sec. HIl but with the electron and hole co- 
ordinates being re-expressed in terms of ri, r2, and the 
fixed vector R. 

The center-of-mass constraint may not be used to cal- 
culate the interaction potential at zero exciton-exciton 
separation, because in that limit the repulsion becomes 
strong enough to dissociate the two individual excitons. 
The ground state of Eq. ([2]) at very small R is thus not 
the physical quantity we require. We have calculated the 
exciton-exciton potential only at separations R for which 
the excitons remain bound. Figure [5] demonstrates this 
effect, exhibiting a potential which decreases at small R 
to physically unreasonable values. 

Our DMC calculations yield a smooth exciton-exciton 
potential. The interaction energies shown in Fig.[5]do not 
deviate from simple dipole-dipole repulsion of the form 



2(f 



(3) 



by more than 3.5 x 10~°R* above an exciton-exciton sep- 
aration of « 7a*g. The fits to the interaction potential 
data are shown in Sec.[Xl The repulsive tails of the inter- 
action (i? > 10a|j) calculated for pairs of excitons with 
0.1 < a < 1 all collapse onto a single curve for each value 
of d when scaled into excitonic units, showing a maximum 
deviation from each other and Eq. ([3]) of 8 x lO^^i?*. 

For an electron/hole mass ratio of cr = 0, our results 
should reduce to the exciton-exciton interaction under 
the heavy-hole approximation4 Figure [S] demonstrates 
the agreement with the interaction potentials calculated 
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FIG. 5: (Color online) Exciton-exciton interaction potential 
Ei{R) as a function of center-of-mass separation R with a = 
1. The solid lines show the fit to the DMC data. Dashed lines 
show the dipole-dipole interaction energy [Eq. (j3])]. 




FIG. 6: (Color online) Exciton-exciton interaction potential 
Ej{R) as a function of hole- hole separation R for the heavy- 
hole case (cr = 0) with d — 0.9, 1.0, and l.la*g. The solid 
lines show the interaction potential from Ref. @. 



by Schindler and Zimmermann. Our points are slightly 
below the curves of Ref. 1^3 at the potential minimum, 
although the agreement is in general very good, and well 
within statistical error at large R. Equation (|A1[) (in 
Sec. El shows a functional form suitable for fitting to 
our data. 

For large layer separations d the interaction is purely 
repulsive, whereas for smaller d the interaction is at- 
tractive at short range. The critical point in the bind- 
ing occurs near the layer separation for which the min- 
imum in the exciton-exciton interaction potential dis- 
appears. Schindler and Zimmermann's approach uses 
a model exciton-exciton interaction potential which de- 
pends on the layer separation d but not the mass ratio cri^ 
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FIG. 7: (Color online) Exciton-exciton interaction potential piQ. 8: (Color online) PDF ^^["'^'''(r) for an isolated electron- 
Ei(R) as a function of (constrained) center-of-mass separa- hole pair from the exact solution of Eq. (2) in Ref.0, shown 
tion 7i for d = 0.9as and ct = and 1. at several layer separations for o" = 1. 



Constraining the center of mass rather than the hole po- 
sitions allows us to observe the interaction potential for 
different mass ratios and layer separations, so that we do 
not need to apply the interaction potential obtained in 
one system to another with different parameters. The 
dependence of the interaction potential upon a is clear 
from Fig. [71 and is consistent with the results shown in 
Fig. m in which bicxcitons are stable at = Q.9a*g for 
cr = but not (7 = 1. 

For the strictly two-dimensional case (rf = 0), we can 
compare our values of the Haynes factor, fn = Eb/Ex, 
with those of previous work. Usukura et al. performed 
numerically exact variational calculations, finding fn = 
0.665 for cr and Jh = 0.193 for a = These 
data agree well with our values of Jh — 0.670(3) and 
0.19287(2) for cr = and 1, respectively. 



IV. PDF'S IN BIEXCITONS 
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FIG. 9: (Color online) Extrapolated electron-electron PDF 
gcT{r) for bound biexcitons with a = 1. The hole-hole and 
electron-electron PDF's are identical for equal electron and 
hole masses. 



The PDF's of electrons and holes in bicxcitons re- 
veal important information about the physics of biexciton 
binding. The electron-electron PDF is defined as 

5cc(r) = ^((S(|reT-rex|-r)), (4) 

where rcf and are the positions of the up- and down- 
spin electrons and the angled brackets denote the average 
over sets of electron and hole coordinates distributed as 
the square of the ground-state wave function. The hole- 
hole PDF is defined in a similar fashion. The electron- 
hole PDF is defined to be 




where roa^ — r{[^^ is the in-plane separation of an elec- 
tron and a hole. The PDF's may be accumulated within 
QMC by binning interparticle distances. The errors in 
the VMC and DMC estimates of the PDF [g^^'^{r) and 
^DMC^^jj g^j.^, lijicar in the error in the trial wave func- 
tion; however, the error in the extrapolated estimate 
(r) = 2gDMC(^-) _ gVMC(^) jg second order in the 

error in the wave function.— Our VMC and DMC PDF's 
are very close to one another, so the errors in our extrap- 
olated estimates are small. The PDF's presented here 
have been normalized such that 

/•OO 

/ 2TTrg''''\r)dr = 1. (6) 
Jo 
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FIG. 10: (Color online) Extrapolated electron-hole PDF 
g'^ir) for the biexciton system with o = 1 and several bilayer 
separations d. 



[TOl the quantity 2glf 



- ath^'^^-r) plotted in Fig. M 
allows one to see the inter-exciton electron-hole PDF su- 
perimposed on the change in the intra-exciton PDF due 
to the presence of the other exciton. The peaks in Fig. [11] 
occur at the same separation as those in Fig. [9l confirm- 
ing that cxcitons retain their identity in bound biexcitons 
for large layer separations, even when electrons and holes 
have equal masses. For zero layer separation there is no 
discernable peak, however, and the function rises sharply 
to a maximum at zero interparticle separation. This may 
be due to the large change in the single-exciton PDF due 
to the presence of the other exciton swamping the inter- 
exciton electron- hole PDF. We are thus unable to con- 
clude with certainty that cxcitons retain their identities 
in bound biexcitons throughout the region of biexciton 
stability. 



V. CONCLUSIONS 
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We have carried out a QMC study of the interaction 
between pairs of cxcitons in bilayer systems. We have cal- 
culated the exciton-exciton interaction potential by con- 
straining the center-of-mass separation, which we believe 
gives a more accurate pair potential at short range than 
the potential calculated by assuming the holes to be in- 
finitely heavy.— We find that for large layer separations, 
cxcitons retain their identity when they bind to form a 
biexciton, suggesting that treating cxcitons as individ- 
ual particles is a reasonable approximation. However, by 
solving the Schrodinger equation for all four particles in 
a biexciton, we find that the range of layer separations 
and mass ratios over which biexcitons are stable is some- 
what larger than the region of stability predicted using 
exciton-exciton pair potentials. 



FIG. 11: (Color online) Biexciton electron-hole PDF relative 
to the single-exciton PDF, 2g°f (r) - gt^^^^'ir), at a = 1 and 
several bilayer separations d. 



Figure [8] shows the electron-hole PDF for a single exci- 
ton, obtained from the exact numerical solution to 
Eq. (2) inRef.0. Figures [51 and [TUl show electron-electron 
and electron- hole PDF's, respectively, for the biexciton 
system with a = 1. At smaller layer separations the 
electron-hole PDF exhibits a larger peak at zero inter- 
particle separation, and decays more rapidly with inter- 
particle distance. 

The size of the biexciton is most easily judged by ex- 
amining the electron-electron PDF (which is identical to 
the hole- hole PDF for cr = 1). The size of the biexciton 
diverges as the critical layer separation (dcrit = 0.43a|j 
for cr = 1) is approached. At zero layer separation, the 
electron-electron PDF is negligible for interparticle dis- 
tances larger than 3a|j and has a maximum at 0.3a|j. 

Although a second peak cannot be discerned in Fig. 
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APPENDIX A: FIT TO THE EXCITON-EXCITON 
POTENTIAL 

The exciton-exciton potential curves with ct = in Fig. 
inimay be fitted to a function of the form, 

/ 1000 »2 ^ P3\ / PiR^\ 

= (,p^ + -?r + i?^ + i?jj°"H"wooj 
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Parameter d = 0.9an d = l.Oao d — Llat 



Pi 


—70.18 


—66.25 


—61.89 


P2 




— ^Ooo 


— ^ou^ 


P3 


4296 


8422 


12560 


Pi 


8.086 


7.319 


6.420 


P5 


21520 


30740 


40610 


P6 


-15100 


-54120 


-98510 


P7 


0.1284 


0.1451 


0.2424 



TABLE I: Coefficients appearing in Eq. (|Aip allowing the re- 
production of fits to the points shown in Fig. |6l Performing 
the fits using data with R > Sag yields errors of 0.79, 1.1 
and 1.4 per data point for d — 0.9, 1.0, and l.la*g, respec- 
tively. 



Parameter d 


= 0.2a^ 


d = 0.5a|3 


d = 0.9a|5 


Pi 


2302 


1463 


-6797 


P2 


-8947 


-12580 


73200 


P3 


1316 


5.813 


24.94 


Pi 


0.1123 


4.703 


2.465 



TABLE IL Coefficients appearing in Eq. (|A2|I allowing the 
reproduction of fits to the DMC results in Fig. [5] Performing 
the fits using data with R > 2.50% yields errors of 1.55, 
and 1.07 per data point for d = 0.2 and 0.5ag, respectively. 
The error is larger for d — 0.90^, the purely repulsive 
curve, but the maximum deviation from the data points is 
only 1.2 X 10"^i?;. 



where d is the layer separation and pi , . . . , P7 are the 
fitting parameters. The function has the correct long- 
range behavior, Ej cx 2(1^ /B? for R oo. The fitting 
parameter values are shown in Table HI 

The interaction potentials in Fig. [5] with tr = 1 may be 
fitted to a function similar to Eq. (jAip . This time the 
form is 



where the long range behavior is once again reproduced 
correctly and each of the terms in the first bracket has 
a physical interpretation. The l/R^ term may be as- 
sociated with quadrupole-quadrupole repulsion and the 
\/R^ term with van der Waals attraction. The signs of 
the fitting parameters are consistent with this interpre- 
tation for d = 0.2 and 0.5a|j. The parameter values are 
shown in Table [TTl 
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